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DISCLAIMER 



Although each program has been tested by its contributor, no 
warranty, express or implied, is made by the contributor or 
1620 USERS Group, as to the accuracy and functioning of the 
program and related program material, nor shall the fact of 
distribution constitute any such warranty, and no responsibility 
is assumed by the contributor or 1620 USERS Group, in con- 
nection therewith. 



1620 USERS GROUP PROGRAM REVIEW AND EVALUATION 
(fill out in typewriter or pencil, do not use ink) 

■ogram No. . Date 

Program Name: 

1. Does the abstract adequately describe what the program is and what Yes No 

it does? 

Comment 

2. Does the program do what the abstract says ? Yes No 

Commen t 

3. Is the Description clear, understandable, and adequate? Yes No 
Comment 

4. Are the Operating Instructions understandable and in sufficient detail? Yes_ No 

Comment 

Are the Sense Switch options adequately described (if applicable)? Yes No 

Are the mnemonic labels identified or sufficiently understandable? Yes No 

Comment 

5. Does the source program compile satisfactorily (if applicable)? Yes No 

Comment 

f\ 

6. Does the object program run satisfactorily? Yes No 

Comment 

7. Number of test cases run . Are any restrictions as to data, 

size, range, etc, covered adequately in description? Yes No 

Comment 

8. Does the Program Meet the minimal standards of the 1620 Users 

Group? Yes No 

Comment ' 

9. Were all necessary parts of the program received? Yes No 

Comment 

10. Please list on the back any suggestions to improve the usefulness of the program. 
These will be passed onto the author for his consideration. 

Please return to: Your Name 



Mr. Richard L. Pratt Company 

Data Corporation 

7500 Old Xenia Pike Address 

j) Dayton, Ohio 45432 

User Group Code 

THIS REVIEW FORM IS PART OF THE 1620 USER GROUP ORGANIZATION'S PROGRAM 
REVIEW AND EVALUATION PROCEDURE. NONMEMBERS ARE CORDIALLY INVITED 
TO PARTICIPATE IN THIS EVALUATION. 
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Computation of Bessel Functions of Integral Order 
by 

C. E. Grosch 

Many problems involving axi- symmetric potentials have 
as a final step the evaluation of infinite series of Bessel 
functions of integral order. That is, the solution may be 
expressed in terms of an infinite series, each term of which 
contains a Bessel function of integral order, or the solu- 
tion is expressible as an integral, the integrand of which 
contains one or more infinite series of Bessel functions. 
Since the summing of the series or the numerical quadrature 
is tedious and time consuming when performed on a desk cal- 
culator there is an obvious need for a computer program to 
generate Bessel functions. This program can then be incor- 
porated in programs to sum the series or do the numerical 
integration. 

Such a program is described in this note. This program 
is based on an analysis presented in Ref. (l). Since Ref. (l) 
contains an error, the corrected analysis will be presented 
below. 

Analysis 

The formula 

.W*) =f J n (x) " J n-l(*) . CD 



- 2 - 



for generating Bessel functions of the first kind for fixed 
x , starting from J D ( X ) and J-^x) is satisfactory so long 
as n<x. If n>x the process becomes unstable. If one 
chooses, however an N sufficiently large compared to x , so 
that j jj+x( x ) = , to some order of approximation and 
J*(x) < < 1 but not zero then the use of the recursion for- 
mula, in a backwards manner is a stable process. 

Since N > x and N > > 1 one can use as a starting 
approximation 2 ^ 

j* (eN) _ e N exp [N(l - e 2 ) V2 ] (a) 

J2ZN(1 - e 2 ) 1 / 4 [1 + JT^f 

< e < 1 
# 

Assume that each J n ( x ) differs from the correct value 
J n (x), i. e. 

J n ( x ) = Vn^ (3) 
Since, to our order of approximation 

<W X > =f V x ) W 

and 

W x ) =f J N< X ) (5) 

*N-1 = *N ^ 

In fact, it is easily shown that 

K = K 1 = =K N-1 =K N = K < 7 > 
Using the identity 

1 = J (x) + 2 . J (x) (8) 

° n=l dn 
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± = J*{x) +2 J (x) (9) 

* ° n=l 2n 

It is now only necessary to determine a minimum N 
large enough so that the process is sufficiently accurate. 
This has been done in Ref . (l) where the following relations 
are given 

10" 51 < x < 10 N — 35.0/(3.5-ln x) (10) 

10 < x < 500 N - (21/20 )x + 25 (ll) 

Range 

Since the power series expansion of J n ( x ) converges 
rapidly for x < < 1 and the asympto'tic series is very accurate 
for x > > 1 , the program was written for the intermediate range 
where neither the power series nor asymptotic series is useful. 
The range chosen for this program Is 0.001 <. x _< 200.0 . 

Accuracy 

The accuracy of this program has been checked by computing 
J n {x) for x = 0.001, 0.100, 1.000, 3.000, 5.000, 7.000, 
10.000, 100.000 and 200.000, and comparing with tabulated 
values. The largest error is +1 in the seventh significant 
digit. 

Running Time 

Approximately five values of J n ( x ) are computed per 
second. Therefore the running time is approximately 0.2N 
sec, that is the running time, T, is, in seconds 



© 
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(7.0/(5.5 
5.0 + 0. 



5-ln x) 0.001 £ x £ 10 
2x 10 £ x <; 200 



Input - Output 

The Input, output for this program is the IBM 1622 
Card Read- Punch unit. 

Likely Changes 

The main use of this program is believed to be as part 
of a larger program which utilizes Bessel Functions in some 
manner. One would then remove the input-output statements. 
Another modification could occur if only J m , J n , Jy ... were 
needed. These can be easily located by noting that 

c(i) = J N+1 

C(2) = J N - 



C(N + 2) = J c 

Then C(N + 2 - m) = . J m and so on. 

Finally, there may be need for the Bessel Functions of 
the second kind Y n (x). These are easily found once the 
J n (x) are known by using^ 1 ^ 

Yjx) = (2/tt) [J n (x) {y + m §) - 2 ^ I J-i) n J 2n (x)] 



r 



2' 

= 0.57721566 



Yl (x) = [J^xjYjx) - ^ ]/J (x) 
and the recurvise relation 

W*> = f V x > " Y n-l< x > 



Description of Program 

The value of x is read in and punched out. Al is equal 
to ^ir) -1 / 2 . The next three statements initialize sum and 
the C 1 . The following four statements, through 4, calculate 
N in floating point and truncate to a fixed point integer. 
The formulas for N , (10 ) and (11 ) are for a minimum N , 
therefore 1.0 has been added to ensure that the minimum is 
exceeded. 



In the program 


C l 


= W x ) 




C 2 


= J N ( X ) 






= J N-l( X > 



C N+2 = J Q (x) 

therefore the final value of i is N+2. Statement 5 computes 

this final value. The next six statements compute C 2 using 

eq. (2). The statements through 6 compute the , i.e. the 

* 

J n (x). The four statements following 6 determine whether or 
not N is odd or even. If N is even then the sum 

N N 
2 C, = 2 J* 
1=2,4,6 1 n=l 2n 

is carried out in statements 7 and 9. The next two statements 

add to this sum 

N+2 N 
2 C< = 2 J* 



If N is odd the corresponding sums (over odd i 3 i.e. even n) 
are calculated in statements 8 through 12. 

In either case, 13 computes K and the statements through 
14, compute n, J n (x) and punch them. The program then cycles 
hack to read another x. 

Sample printouts are appended. 
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C COMPUTATION OF BESSEL FUNCT I ONS , F I RST KIND* INTEGRAL ORDER 

C FOR ARGUEMENTS IN THE RANGE ♦ GREATER THAN • 00 1 « LESS THAN 200* 

99 FORMAT (F13.3) 
98 FORMAT U9.E16.8) 

DIMENSION C(240> 
14 READ 99.X 
PUNCH 99.X 
Al =0.39894240 
SUM=0.0 
DO 2 1=1 .240 

2 Cf I ) =0.0 
IF (X-10.0) 3.3.4 

3 M=<3R.0/< 3.5-LOGCX) ) )+l .O 
GO TO 5 

4 N=l .0fi*X+?6.0 

5 IFNO=N+2 
FLN = N 
F=X/FLN 
Y=l .0-E*F 
A=SORT(Y) 

A2=< F**FLN>*< (FLN*A) **( -0.5 ) )*( ( 1 ,0 + A)** (-FLN) ) 
C< 2) =A1*A2*EXP(FLN*A ) 
Xl=2.0/X 
DO 6 I = 1 . N 
M=N-I+1 
FLM=M 

6 C ( I +2 ) =FLM*x 1 *C < I + 1 > -C < I ) 
D=O.B*FLN 

J=D 
FLJ = J 

IF <D-FLJ) 7.7,8 

7 K = 2 
GO TO 9 

8 K=l 

9 DO 10 l=K,N.2 

10 SUM=SUM+C( I ) 
DO 11 I=K,IEND.2 

11 SUM=SUM+C(I) 
FLK=1.0/SUM 
DO 12 1=1 , I END 

12 C( I )=FLK*C( I ) 
DO 13 I =1.1 END 
NP=N-I+2 

13 PUNCH 98.NP,C(I) 
GO -TO 14 
END 
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